function handles = analyzeArrayOfAveragedTraces(handles)
% analyzeArrayOfAveragedTraces
%
% gs june 2005 
% ------------------------------------------------------------------------------

% ANALYSIS METHOD
% analysisMethod = get(handles.methodInputMap, 'Value');
analysisMethod = 3;
% method 1 -- old way: baseline-subtract using MEAN; divide baseline into baseline/response
% method 2 -- new way: baseline-subtract using MEDIAN; simply calc SD of baselines
% method 3 -- variant: baseline-subtract using MODE; calc SD of baselines

% DATA etc
dataArray = handles.data.map.mapActive.bsArray; % dataArray
[rows, numTraces] = size(dataArray); % must be 2D data array; col num = trace num
sr = handles.data.acq.sampleRate;

mapPattern = handles.data.map.mapActive.mapPattern;
traceData = []; % unused?

% % CELL parameters
% if get(handles.chkWithRsRm, 'Value')
%     if ~handles.data.analysis.currentClamp
%         handles.data.analysis.skipCm = 0; % option here to skip Cm calc, e.g. if it's too slow
%         handles = arrayCellParams(handles, dataArray);
%     else
%         disp('Cell parameter calculations not yet implemented for current clamp.');
%     end
% end

% TIME vector
totalTime = (rows-1)/sr; % in sec, not msec
xTimeAxis = linspace(0, totalTime, rows)';
handles.data.map.mapActive.xTimeAxis = xTimeAxis;
minX = 0; % unused?
maxX = max(xTimeAxis); % unused?

% BASELINE settings, mean, sd, baseline-subtracted array
if analysisMethod == 1
    baselineStartIndex = round(handles.data.analysis.baselineStart*sr)+1;
    baselineEndIndex = round(handles.data.analysis.baselineEnd*sr)+1;
    baselineMeans = mean(dataArray(baselineStartIndex:baselineEndIndex,:));
    handles.data.map.mapActive.baselineMeans = baselineMeans';
    baselineSDs = std(dataArray(baselineStartIndex:baselineEndIndex,:));
    handles.data.map.mapActive.baselineSDs = baselineSDs';
    % % baseline-subtracted array:
    % baselineArray = repmat(baselineMeans, rows, 1);
    % bsArray = dataArray - baselineArray; 
    % handles.data.map.mapActive.bsArray = bsArray;
    bsArray = handles.data.map.mapActive.bsArray;
    
    % NOISE of baseline
    %       b0 = 0-50 msec, baseline for baseline window
    %       r0 = 50-100 msec, baseline window
    %       b1 = 50-100 msec, baseline for synaptic window
    %       r1 = 110-160 msec, synaptic window
    % % mean 
    % b0 = mean(state.uncaging.map.bsArray(1:500, :));
    % b1 = mean(state.uncaging.map.bsArray(500:1000, :));
    % r0 = mean(state.uncaging.map.bsArray(500:1000, :));
    % r1 = mean(state.uncaging.map.bsArray(1101:1601, :));
    % peak 
    b0 = min(handles.data.map.mapActive.bsArray(1:500, :));
    b1 = min(handles.data.map.mapActive.bsArray(500:1000, :));
    r0 = min(handles.data.map.mapActive.bsArray(500:1000, :));
    r1 = min(handles.data.map.mapActive.bsArray(1101:1601, :));
    % 
    handles.data.map.mapActive.ampHistoBaselineData = -(r0-b0);
    handles.data.map.mapActive.histoBaseSDs = std(handles.data.map.mapActive.ampHistoBaselineData);
    handles.data.map.mapActive.ampHistoSynapticData = -(r1-b1);
    
elseif analysisMethod == 2
    % obtain median values of each baseline region
    baselineStartIndex = round(handles.data.analysis.baselineStart*sr)+1;
    baselineEndIndex = round(handles.data.analysis.baselineEnd*sr)+1;
    baselineMedians = median(dataArray(baselineStartIndex:baselineEndIndex,:));
    handles.data.map.mapActive.baselineMedians = baselineMedians';
    
    % obtain the SD in the baseline region
    baselineSDs = std(dataArray(baselineStartIndex:baselineEndIndex,:));
    handles.data.map.mapActive.baselineSDs = baselineSDs';
    
    % baseline-subtracted array:
    baselineArray = repmat(baselineMedians, rows, 1);
    bsArray = dataArray - baselineArray; 
    handles.data.map.mapActive.bsArray = bsArray;
        
elseif analysisMethod == 3 
    % obtain mode values of each baseline region
    handles.data.map.mapActive.baselineModes = [];
    for i = 1 : numTraces
        baselineStartIndex = round(handles.data.analysis.baselineStart*sr)+1;
        baselineEndIndex = round(handles.data.analysis.baselineEnd*sr)+1;
        data = dataArray(baselineStartIndex:baselineEndIndex, i);
        handles.data.map.mapActive.baselineModes(i) = mode(data);
    end
    
    % obtain the SD in the baseline region
    baselineSDs = std(dataArray(baselineStartIndex:baselineEndIndex,:));
    handles.data.map.mapActive.baselineSDs = baselineSDs';
    
    % baseline-subtracted array:
    baselineArray = repmat(handles.data.map.mapActive.baselineModes, rows, 1);
    bsArray = dataArray - baselineArray; 
    handles.data.map.mapActive.bsArray = bsArray;
end

if analysisMethod == 2 | analysisMethod == 3
        % NOISE of baseline
    %       b0 = 0-50 msec, baseline for baseline window
    %       r0 = 50-100 msec, baseline window
    %       b1 = 50-100 msec, baseline for synaptic window
    %       r1 = 110-160 msec, synaptic window
    % % mean 
    % b0 = mean(state.uncaging.map.bsArray(1:500, :));
    % b1 = mean(state.uncaging.map.bsArray(500:1000, :));
    % r0 = mean(state.uncaging.map.bsArray(500:1000, :));
    % r1 = mean(state.uncaging.map.bsArray(1101:1601, :));
    % peak 
    b0 = min(handles.data.map.mapActive.bsArray(1:500, :));      % these are all vectors
    b1 = min(handles.data.map.mapActive.bsArray(500:1000, :));
    r0 = min(handles.data.map.mapActive.bsArray(500:1000, :));
    r1 = min(handles.data.map.mapActive.bsArray(1101:1601, :));
    % 
    handles.data.map.mapActive.ampHistoBaselineData = -(r0-b0);  % vector
    handles.data.map.mapActive.histoBaseSDs = std(handles.data.map.mapActive.ampHistoBaselineData); % single value (not a vector)
    
    % new way
    handles.data.map.mapActive.ampHistoBaselineData = -(r0-b0);
    
    handles.data.map.mapActive.histoBaseSDs = mean(std(dataArray(baselineStartIndex:baselineEndIndex,:)));
    handles.data.map.mapActive.ampHistoSynapticData = -(r1-b1);
end

% STIMULUS timing
stimOn = handles.data.analysis.stimOn;
stimOnInd = round(stimOn*sr+1);
handles.data.analysis.stimOnInd = stimOnInd;

% WINDOWS *******************************************************************
directWindowStart = handles.data.analysis.directWindowStart;
directStartInd = round(directWindowStart * sr + 1);
directWindow = handles.data.analysis.directWindow;
directEndInd = directStartInd + round(directWindow * sr);
synStart = directWindow;
synStartInd = directEndInd;
synapticWindow = handles.data.analysis.synapticWindow;
synEndInd = directStartInd + round(synapticWindow * sr);

% THRESHOLDS ****************************************************************
% NB: these are vectors
synLevel = handles.data.analysis.synThreshold; 
dirLevel = synLevel; 
if analysisMethod == 1 % method 1 -- based on baseline mean
    handles.data.map.mapActive.dirPosThresh = handles.data.map.mapActive.baselineMeans + dirLevel * handles.data.map.mapActive.histoBaseSDs;
    handles.data.map.mapActive.dirNegThresh = handles.data.map.mapActive.baselineMeans - dirLevel * handles.data.map.mapActive.histoBaseSDs;
    handles.data.map.mapActive.synPosThresh = handles.data.map.mapActive.baselineMeans + synLevel * handles.data.map.mapActive.histoBaseSDs;
    handles.data.map.mapActive.synNegThresh = handles.data.map.mapActive.baselineMeans - synLevel * handles.data.map.mapActive.histoBaseSDs;
elseif analysisMethod == 2 | analysisMethod == 3 % 
    handles.data.map.mapActive.dirPosThresh = dirLevel * handles.data.map.mapActive.histoBaseSDs;
    handles.data.map.mapActive.dirNegThresh = -dirLevel * handles.data.map.mapActive.histoBaseSDs;
    handles.data.map.mapActive.synPosThresh = synLevel * handles.data.map.mapActive.histoBaseSDs;
    handles.data.map.mapActive.synNegThresh = -synLevel * handles.data.map.mapActive.histoBaseSDs;
end

handles.data.map.mapActive.synThreshpAmV = synLevel * handles.data.map.mapActive.histoBaseSDs;
handles.data.map.mapActive.baselineSD = handles.data.map.mapActive.histoBaseSDs;

% RESPONSES *****************************************************************
% strategy: measure response params over entire post-stim window (direct + syn)
% aftwerward, determine whether responses were direct or syn

% MEAN
genMeanAmp = mean(bsArray(stimOnInd:synEndInd,:));
handles.data.map.mapActive.genMeanAmp = genMeanAmp';

% MIN
[genMinAmp, genMinPkLatInd] = min(bsArray(stimOnInd:synEndInd,:));
genLatOfMin = (genMinPkLatInd-1)/sr;
handles.data.map.mapActive.genMinAmp = genMinAmp';
handles.data.map.mapActive.genLatOfMin = genLatOfMin';

% MAX
[genMaxAmp, genMaxPkLatInd] = max(bsArray(stimOnInd:synEndInd,:));
genLatOfMax = (genMaxPkLatInd-1)/sr;
handles.data.map.mapActive.genMaxAmp = genMaxAmp';
handles.data.map.mapActive.genLatOfMax = genLatOfMax';

% INTEGRAL
genIntegral = trapz(bsArray(stimOnInd:synEndInd,:))/sr;
handles.data.map.mapActive.genIntegral = genIntegral';

% THRESHOLDED MIN
result = arrayEventFinder(dataArray(stimOnInd:synEndInd,:), handles.data.map.mapActive.synNegThresh, 'negative');
events = zeros(numTraces,1);
if ~isempty(result)
    tracesWithEvents = result(:,1);
    timeIndicesOfEvents = result(:,2);
    onsetLatEvents = result(:,2)/sr;
    events(tracesWithEvents) = onsetLatEvents;
end
handles.data.map.mapActive.genLatOfMinOnset = events;

% THRESHOLDED MIN -- spontaneous event analysis of baseline interval
% TODO: finish this; event finder must be able to count number of events
result = arrayEventFinder(dataArray(baselineStartIndex:baselineEndIndex,:), handles.data.map.mapActive.synNegThresh, 'negative');
events = zeros(numTraces,1);
if ~isempty(result)
    tracesWithEvents = result(:,1);
%     timeIndicesOfEvents = result(:,2);
%     onsetLatEvents = result(:,2)/sr;
%     events(tracesWithEvents) = onsetLatEvents;
end
handles.data.map.mapActive.spontEventNum = numel(tracesWithEvents);
totalBaselineTime = numTraces * (handles.data.analysis.baselineEnd - handles.data.analysis.baselineStart);
handles.data.map.mapActive.spontEventRate = numel(tracesWithEvents) / totalBaselineTime;

% THRESHOLDED MAX
result = arrayEventFinder(dataArray(stimOnInd:synEndInd,:), handles.data.map.mapActive.synPosThresh, 'positive');
events = zeros(numTraces,1);
if ~isempty(result)
    tracesWithEvents = result(:,1);
    timeIndicesOfEvents = result(:,2);
    onsetLatEvents = result(:,2)/sr;
    events(tracesWithEvents) = onsetLatEvents;
end
handles.data.map.mapActive.genLatOfMaxOnset = events;

% % % % TYPES OF RESPONSES according to clamp mode, response polarity
% % % mode = state.uncaging.analysis.currentClamp;
% % % switch mode
% % % case 0 % VC
    % assumes typical negative Vh -- inward currents are excitatory, outward currents inhibitory

    % first the direct responses: suprathreshold excitatory responses starting within direct window
    handles.data.map.mapActive.edBinary = handles.data.map.mapActive.genLatOfMinOnset > (directWindowStart - stimOn) & ...
        handles.data.map.mapActive.genLatOfMinOnset <= directWindow;
    handles.data.map.mapActive.edAmp = handles.data.map.mapActive.genMinAmp.*handles.data.map.mapActive.edBinary;
    handles.data.map.mapActive.edDelay = handles.data.map.mapActive.genLatOfMin.*handles.data.map.mapActive.edBinary;
    handles.data.map.mapActive.edDelayOnset = handles.data.map.mapActive.genLatOfMinOnset.*handles.data.map.mapActive.edBinary;
    
    % now the synaptic responses -- suprathreshold excitatory responses starting within synaptic window
    handles.data.map.mapActive.ecBinary = handles.data.map.mapActive.genLatOfMinOnset > synStart & ...
        handles.data.map.mapActive.genLatOfMinOnset <= synapticWindow;
    handles.data.map.mapActive.ecAmp = handles.data.map.mapActive.genMinAmp.*handles.data.map.mapActive.ecBinary;
    handles.data.map.mapActive.ecDelay = handles.data.map.mapActive.genLatOfMin.*handles.data.map.mapActive.ecBinary;
    handles.data.map.mapActive.ecDelayOnset = handles.data.map.mapActive.genLatOfMinOnset.*handles.data.map.mapActive.ecBinary;
    
    handles.data.map.mapActive.icBinary = handles.data.map.mapActive.genLatOfMaxOnset > synStart & ...
        handles.data.map.mapActive.genLatOfMaxOnset <= synapticWindow;
    handles.data.map.mapActive.icAmp = handles.data.map.mapActive.genMaxAmp.*handles.data.map.mapActive.icBinary;
    handles.data.map.mapActive.icDelay = handles.data.map.mapActive.genLatOfMax.*handles.data.map.mapActive.icBinary;
    handles.data.map.mapActive.icDelayOnset = handles.data.map.mapActive.genLatOfMaxOnset.*handles.data.map.mapActive.icBinary;
    
    %integrals
    
% % % case 1 % CC
% % %     % assumes typical negative Vr -- positive responses excitatory, negative responses inhibitory
% % % 
% % %     % first the direct responses: suprathreshold excitatory responses starting within direct window
% % %     state.uncaging.map.edBinary = state.uncaging.map.genLatOfMaxOnset > (directWindowStart - stimOn) & ...
% % %         state.uncaging.map.genLatOfMaxOnset <= directWindow;
% % %     state.uncaging.map.edAmp = state.uncaging.map.genMaxAmp.*state.uncaging.map.edBinary;
% % %     state.uncaging.map.edDelay = state.uncaging.map.genLatOfMax.*state.uncaging.map.edBinary;
% % %     state.uncaging.map.edDelayOnset = state.uncaging.map.genLatOfMaxOnset.*state.uncaging.map.edBinary;
% % %     
% % %     % now the synaptic responses -- suprathreshold excitatory responses starting within synaptic window
% % %     state.uncaging.map.epBinary = state.uncaging.map.genLatOfMaxOnset > synStart & ...
% % %         state.uncaging.map.genLatOfMaxOnset <= synapticWindow;
% % %     state.uncaging.map.epAmp = state.uncaging.map.genMaxAmp.*state.uncaging.map.epBinary;
% % %     state.uncaging.map.epDelay = state.uncaging.map.genLatOfMax.*state.uncaging.map.epBinary;
% % %     state.uncaging.map.epDelayOnset = state.uncaging.map.genLatOfMaxOnset.*state.uncaging.map.epBinary;
% % %     
% % %     state.uncaging.map.ipBinary = state.uncaging.map.genLatOfMinOnset > synStart & ...
% % %         state.uncaging.map.genLatOfMinOnset <= synapticWindow;
% % %     state.uncaging.map.ipAmp = state.uncaging.map.genMinAmp.*state.uncaging.map.ipBinary;
% % %     state.uncaging.map.ipDelay = state.uncaging.map.genLatOfMin.*state.uncaging.map.ipBinary;
% % %     state.uncaging.map.ipDelayOnset = state.uncaging.map.genLatOfMinOnset.*state.uncaging.map.ipBinary;
% % %     
% % % case 2 
% % % end
                
% assign responseType based on binary vectors for ed, ec/ep, ic/ip, nr
% % % if ~state.uncaging.analysis.currentClamp
    bitArray = [num2str(handles.data.map.mapActive.icBinary) ...
            num2str(handles.data.map.mapActive.ecBinary) ...
            num2str(handles.data.map.mapActive.edBinary)];
    codes = bin2dec(bitArray);
    translator = ['NR'; 'ED'; 'EC'; '  '; 'IC'; '  '; 'EC'; '  '];
    handles.data.map.mapActive.responseType = translator(codes+1,:);
% % % else
% % %     bitArray = [num2str(state.uncaging.map.ipBinary) ...
% % %             num2str(state.uncaging.map.epBinary) ...
% % %             num2str(state.uncaging.map.apBinary)];
% % %     codes = bin2dec(bitArray);
% % %     translator = ['NR'; 'AP'; 'EP'; '  '; 'IP'; '  '; 'EP'; '  '];
% % %     state.uncaging.map.responseType = translator(codes+1,:);
% % % end

handles = makeMapDataArray(handles);
handles = array2Dplot(handles);



% ============ subfunctions ===================
function val = mode(data);
sorted = sort(data);
[dummy, idx] = unique(sorted);
[dummy, idx] = max(diff(idx));
val = sorted(idx);
